#rm(list=ls())
library(plyr)
library(ggplot2)

#### Replication Code for Grossman, Phillips and Rosenzeig (2017), "Opportunistic Accountability: State-Society Bargaining Over Shared Interests" 

#### As agreed with the editors of Comparative Political Studies, one of the two datasets, (file "UNICEF_compliance_data_2012_to_2015_processed"), is not being currently shared as that data is owned by a third party.

#setwd("~/Dropbox/Polio Paper (1)/Content for CPS paper/code/Replicable Code and Figures")
#setwd("C:/Users/Jonny/Dropbox/Projects/Polio Paper/Content for CPS paper/code/Replicable Code and Figures")

#Load Data
cases <- read.csv("Data/Polio Incidence Data/UNICEF_polio_incidence_nigeria_2000_2013.csv", header=TRUE, stringsAsFactors=FALSE)
data <- read.csv("Data/Compliance Data/UNICEF_compliance_data_2012_to_2015_processed.csv", check.names=FALSE)

#### Figure 1: Polio cases in Nigeria have dropped since 2009 ####
cases$DateOfOnset <- as.Date(cases$DateOfOnset,format="%d/%m/%Y")

pdf("Figures/Incidence3.pdf")
hist(cases$DateOfOnset, "months",format="%m/%y",main="Distribution of Polio Cases",xlab="Date (Month/Yr)",ylab="Frequency",col="light blue",freq=T,axes=T)
abline(v=as.Date("0009-06-01"))
abline(v=as.Date("0003-10-01"))
text(as.Date("0002-06-01"),45, "3 northern \ngovernors ban \nvaccine")
text(as.Date("0010-11-01"),35, "Traditional \nrulers brought \non board")
dev.off()

#### Figure 2: Proportion of households refusing to be vaccinated are less than 1% since the end of 2012 ####
Reasons <- c("Household not in microplan", "Household in microplan but not visited", "Child at playground", "Child at social event", "Child at market", "Child at Farm", "Child at school", "Child sick", "Too many rounds", "Fear of OPV side effects", "Polio is rare", "Polio has cure", "There are other remedies available", "Religious belief", "Political differences", "Unhappy with vaccination team", "No pluses given", "No felt need", "No caregiver consent", "Security")
Reasons_Explicit_NC <- c("Too many rounds", "Fear of OPV side effects", "Polio is rare", "Polio has cure", "There are other remedies available", "Religious belief", "Political differences", "Unhappy with vaccination team", "No pluses given", "No felt need", "No caregiver consent")

#Explicit Non-compliance Rate
data$ExplicitNCRate <- rowSums(data[,Reasons_Explicit_NC],na.rm=T)/(data$Marked0to59+data$UnImmun0to59)*100
data$DateMonitor <- as.Date(data$DateMonitor,format="%m/%d/%Y")

data_month <- ddply(data, ~DateMonitor, summarise, ExplicitNCRate=mean(ExplicitNCRate,na.rm=T))

pdf(file="Figures/compliance_exit_voice.pdf", height=5, width=14)
plot(data_month$DateMonitor,data_month$ExplicitNCRate,type="l",xlab="Date (Month-Year)",ylab="Proportion of Households (%)",main="",ylim=c(0,0.8),xaxt='n',lwd=3)
lines(data_month$DateMonitor, data_month$ExplicitNCRate, type="l", col="dark green",lwd=3)
axis(side=1,at=data_month$DateMonitor, labels=format(data_month$DateMonitor, "%m-%Y"))
dev.off()

#### Figure 3: Self-reported Reasons for Household-Level Non-Compliance, January 2013 - September 2015 ####
Reasons_NC <- c("Child at playground", "Child at social event", "Child at market", "Child at Farm", "Child at school", "Child sick", "Too many rounds", "Fear of OPV side effects", "Polio is rare", "Polio has cure", "There are other remedies available", "Religious belief", "Political differences", "Unhappy with vaccination team", "No pluses given", "No felt need", "No caregiver consent")

reasons.d <- as.data.frame(colSums(data[,Reasons_NC],na.rm=T))
names(reasons.d) <- "frequency"
reasons.d$reason <- rownames(reasons.d)
reasons.d$Type<-c(rep("Potentially Hidden Non-Compliance",6),"Non-Compliance Due to Intensity of Polio Campaign",rep("Non-Compliance",8),"Non-Compliance Due to Intensity of Polio Campaign","Potentially Hidden Non-Compliance")
reasons.d$reason <- factor(reasons.d$reason, levels=c("Child at playground","No caregiver consent","Child at Farm","Child at school","Child at social event","Child at market","Child sick","Polio has cure","Political differences","Religious belief","Polio is rare","Fear of OPV side effects","There are other remedies available","No pluses given","Unhappy with vaccination team","Too many rounds","No felt need"))
reasons.d$Type <- factor(reasons.d$Type, c("Non-Compliance Due to Intensity of Polio Campaign","Non-Compliance", "Potentially Hidden Non-Compliance"))

reasons.d <- reasons.d[!(reasons.d$Type %in% "Potentially Hidden Non-Compliance"),]
reasons.d$Type2 <- c("Compliance NOT Costly","Compliance Costly","Compliance Costly","Compliance Costly","Compliance Costly","Compliance Costly","Compliance NOT Costly","Compliance Costly","Compliance NOT Costly","Compliance NOT Costly")
reasons.d$reason <- factor(reasons.d$reason, levels=c("Child at playground","No caregiver consent","Child at Farm","Child at school","Child at social event","Child at market","Child sick","Polio has cure","Religious belief","Polio is rare","Fear of OPV side effects","There are other remedies available","Unhappy with vaccination team","Political differences","No pluses given","Too many rounds","No felt need"))

reasons.d$Type2 <- factor(reasons.d$Type2, levels=c("Compliance NOT Costly","Compliance Costly"))

pdf(file="Figures/Reasons for non-compliance, corresponds to Intervention description and justification section3.pdf", height=5, width=14)
mar.default <- c(5,4,4,5) + 0.1
par(mar = mar.default + c(0, 8, -4, 8)) 
ggplot(reasons.d, aes(x=reason, y= frequency, fill=Type2))+
  geom_bar(stat="identity", position=position_dodge())+
  coord_flip()+
  labs(y="Frequency", x="")+
  scale_y_continuous(expand=c(.001,0),limits=c(0,max(reasons.d$frequency)*1.05),breaks=seq(0, 40000, 5000)) + 
  scale_fill_grey() +
  theme(axis.title.y=element_blank())+
  # ggtitle("Reasons Given by Households for Non-Compliance (9/2014-9/2015)")  + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        legend.text=element_text(size=15),
        text = element_text(size=15),
        legend.title=element_blank())
dev.off()

#### Figure 4: Non-compliance is highest in the most rural and the most urban areas ####
data$NCcount<-rowSums(data[,Reasons],na.rm=T)
data$NCrate<-data$NCcount/(data$Marked0to59+data$UnImmun0to59)*100

data$africa061_log <- log(data$africa061)
data$africa061_log[which(is.infinite(data$africa061_log))] <- NA

ggp.urb<-ggplot(data=data, aes(x=africa061_log, y=NCrate)) + 
  stat_smooth(method="lm", fullrange=F,colour="grey",fill="grey",alpha=.2)+
  stat_smooth(method="lm", fullrange=F,colour="black",fill="black",alpha=.2, formula = y ~ x + I(x^2))+
  geom_point(size=0.1, alpha=0.2) +
  coord_cartesian(ylim=c(0,5)) +
  xlab("Urbanization") + 
  ylab("% Non-compliance") +
  theme_classic()
ggp.urb
ggsave("Figures/Non-Compliance and Urbanization Relationship.png",dpi=600)
